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The history of the carbon cycle is punctuated by enigmatic tran¬ 
sient changes in the ocean's store of carbon. Mass extinction 
is always accompanied by such a disruption, but most disrup¬ 
tions are relatively benign. The less calamitous group exhibits 
a characteristic rate of change whereas greater surges accom¬ 
pany mass extinctions. To better understand these observations, 
I formulate and analyze a mathematical model that suggests 
that disruptions are initiated by perturbation of a permanently 
stable steady state beyond a threshold. The ensuing excita¬ 
tion exhibits the characteristic surge of real disruptions. In this 
view, the magnitude and timescale of the disruption are prop¬ 
erties of the carbon cycle itself rather than its perturbation. 
Surges associated with mass extinction, however, require addi¬ 
tional inputs from external sources such as massive volcanism. 
Surges are excited when CO 2 enters the oceans at a flux that 
exceeds a threshold. The threshold depends on the duration of 
the injection. For injections lasting a time t,-> 10, 000 y in the 
modern carbon cycle, the threshold flux is constant; for smaller 
tj, the threshold scales like tr 1 . Consequently the unusually 
strong but geologically brief duration of modern anthropogenic 
oceanic CO 2 uptake is roughly equivalent, in terms of its potential 
to excite a major disruption, to relatively weak but longer- 
lived perturbations associated with massive volcanism in the 
geologic past. 

carbon cycle | mass extinctions | excitable systems | 
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Introduction 

Earth’s carbon cycle is a loop between photosynthesis, which 
converts carbon dioxide (CO 2 ) to organic carbon, and respira¬ 
tion, which converts organic carbon back to CO 2 (1-4). The 
cycle has undergone many disruptions throughout Earth’s his¬ 
tory (5-7). These events are expressed in the geologic record as 
relatively abrupt and large changes in the isotopic composition of 
sedimentary carbon (8) compared with background fluctuations. 
Fig. 1 shows 2 examples in one 800-ky time series. Such events 
are typically attributed to changes in the fluxes and concentra¬ 
tions of carbon (8-10). Thus, Fig. 1 could reflect the transient 
addition of 2 pulses of isotopically light carbon originating, for 
example, from the respiration of a previously unreactive reser¬ 
voir of organic carbon (10, 11). Such disruptions have also been 
interpreted as geochemical responses to other sources of envi¬ 
ronmental change, including variations in Earth’s orbital param¬ 
eters (12); dissociation of methane hydrate (13); bolide impacts 
(14); biogeochemical innovations (15, 16); and changes in chem¬ 
ical weathering (17), organic carbon burial (18), and volcanic 
emissions (19). These interpretations typically treat the marine 
carbon cycle as a passive recorder of an externally imposed 
stress. The response imprinted in the geochemical record is then 
implicitly assumed to be proportional to the magnitude of the 
forcing (8, 16, 20). 

Given the variety of possible forcing mechanisms, it comes 
as no surprise that the size and timescale of disruptions vary 
immensely, by 2 orders of magnitude during the Phanerozoic 
(0 to 542 Ma) (7). Yet the same data show that major fluctua¬ 
tions of the carbon cycle exhibit a characteristic rate of change. 
Events associated with mass extinction tend to exceed the char¬ 


acteristic rate, whereas others appear bounded by it (7). It seems 
unlikely that a rich diversity of stressors, expressed as propor¬ 
tionate responses in the geochemical record, would exhibit such 
uniformity. 

Here I formulate and analyze an elementary mathematical 
model that shows how the characteristic rate can instead emerge 
within the carbon cycle. The model portrays the upper ocean as 
a chemical reactor open to an incoming flux of dissolved cal¬ 
cium carbonate from rivers and an outgoing flux representing 
carbonate burial in sediments. Within the reactor, the concen¬ 
trations of carbonate species respond not only to imbalances 
in inputs and outputs, but also to imbalances in the biological 
consumption and production of CO 2 . This framework reveals a 
mechanism for autocatalytic amplification of a small but finite 
perturbation of a stable steady state followed by relaxation back 
to the same steady state. The process is analogous to the exci¬ 
tation of an action potential (nerve impulse) in a neuron (21). 
Here, excitations manifest themselves as transient ocean acidifi¬ 
cation; i.e., a temporary increase in the concentration of carbon 
dioxide in the upper ocean. This change is distinguished by a 
characteristic rate. 

These results suggest that the magnitude of a carbon cycle’s 
disruption is determined not by the strength of the cycle’s per¬ 
turbation but rather by the intrinsic dynamics of the system itself. 
Once the addition of CO 2 to the oceans passes a threshold, the 
rate of amplification and the eventual severity of the resulting 
environmental change should be independent of the detailed his¬ 
tory of the perturbation. Moreover, the model indicates that the 
consequences of fast forcing at human timescales may be sim¬ 
ilar to the outcome of slow forcing at geologic timescales (7). 
Within this framework, the greater rates of change associated 
with mass extinction events are straightforwardly interpreted as 
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Fig. 1. Fluctuations of the isotopic composition of carbonate carbon (5 ,3 C) 
during the Eocene period, about 54 Ma (26). Time advances to the right 
and is given with respect to the minimum of the first abrupt downswing, 
an event known as Eocene Thermal Maximum 2 (ETM2) or HI. The sec¬ 
ond event, about 100 ky later, is called H2. The timescale is derived from 
astrochronology (26). 


the consequence of forcing sustained beyond the threshold, due, 
for example, to massive volcanism (22-25). 

This paper is organized as follows. The first section presents 
a 2 -component dynamical system that represents important fea¬ 
tures of the marine carbonate system. Next, this paper analyzes 
the stability of the model’s steady state. Three dynamical regimes 
receive focus: The limit cycle that results when the steady state 
is unstable, the excitations that result when the steady state is 
stable, and the model’s behavior when both the steady state 
and the limit cycle are stable. Particular attention is given to 
the case of excitations and their size and timescale. This paper 
then discusses the implications of these findings for interpret¬ 
ing the geochemical record of past disruptions and predicting 
the carbon cycle’s response to modern anthropogenic perturba¬ 
tions. The conclusion highlights strengths and weaknesses of this 
paper’s principal contributions and points the way toward future 
progress. 

Two-Component Dynamical System 

This section constructs a mathematical model to help explore 
how instabilities might occur in the marine carbon cycle. 
Among the many potentially relevant feedback mechanisms 
(27), I focus on the interactions of shallow-ocean respira¬ 
tion with fluxes of carbon into and out of the shallow ocean. 
These feedbacks are expressed below as a dynamical system 
composed of 2 ordinary differential equations. The objec¬ 
tive is to provide a rational context within which one can 
show how mechanisms operating within the carbon cycle can 
give rise to general features of its past disruptions. Such 
understanding follows from the identification of the ways in 
which the qualitative dynamics of the system—e.g., its sta¬ 
bility and bifurcations—can derive from its intrinsic interac¬ 
tions. The theory of dynamical systems makes such connections 
possible (28-32). 

Although these objectives require the identification of relevant 
mechanisms, they do not necessitate detailed parameterizations. 
The developments below nevertheless give significant attention 
to the way in which specific processes operate in the real carbon 
cycle. I begin by stating essential aspects of carbonate chem¬ 
istry. The model is then stated in its most general form, after 
which pertinent mechanisms are expressed mathematically. A 
final change of variables then leads to a model suitable for 
classical stability analysis. 

The Carbonate System. The model expresses the evolution of 
total alkalinity and the total dissolved inorganic carbon. Total 
alkalinity is defined by the sum (2, 33) 


o = [HCO 3 ] + 2[CO§"] + [B(OH)-] + [OH-] - [H+], [1] 

Total alkalinity represents the excess of bases over acids. Its 
dominant contributions come from bicarbonate ions (HC07) 
and carbonate ions (CO 3 - ); the factor of 2 attached to the 
CO 3 - concentration derives from its double-negative charge 
(33). Total dissolved inorganic carbon (DIC) is symbolized by 
w, mnemonically recalled as the “whole” of inorganic carbon: 

W = [C0 2 ] + [HC0 3 -] + [C0^]. [2] 

The carbonate system is taken to always be in thermal equilib¬ 
rium according to the reactions (2, 33) 

CO 2 + H 2 O —HC 0 J+H+^C 03 ~+ 2 H + . [3] 

The 2 conserved quantities (a and w ) and the 2 reactions of Eq. 
3 specify the equilibrium of the carbonate system. There are 6 
unknowns: Not only a and w, but also the concentrations of CO 2 , 
CO 3 - , HC07, and H + . Specification of any 2 of these quantities 
therefore provides the other 4 via equations (1-3, 33). 

General Formulation. I consider a well-mixed open system in 
which dissolved calcium carbonate (CaC0 3 ) is delivered to the 
shallow ocean by rivers and output from the shallow ocean ulti¬ 
mately into sediments. The assumption of a well-mixed system 
implies that the model addresses timescales longer than the 
timescale of global ocean circulation [about 1,000 y (2, 3)]. The 
precise thickness of the shallow layer is unimportant; however, it 
must be small compared with the average depth of the oceans. 
This assumption allows exchanges of carbon between the shal¬ 
low ocean and the deep ocean and sediments to be expressed in 
analogy to the exchange of heat between a small system and a 
large heat bath. In other words, the concentration of carbon fluc¬ 
tuates within the shallow layer but the mass of carbon below it is 
effectively constant. Thus, the model explicitly considers only the 
chemical state of the upper ocean. 

Fig. 2 illustrates the model. Its mathematical form reads 

a = 2[j in - B(a,w)] [4] 

w = (1 + i/)j in - B(a, w ) + R(a, w) - {w - w 0 )/r w , [5] 

where dots represent time derivatives. The concentrations a and 
w are expressed in units of pmol-kg^ 1 . The factor of 2 arises 
because 2 molar equivalents of alkalinity are associated with each 
mole of CaC0 3 (33). Expressing the model in terms of the con¬ 
served quantities a and w rather than, e.g., the concentrations 
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Fig. 2. Schematic diagram of the model (not drawn to scale). Left and Right 
panels represent, respectively, the evolution of total alkalinity, a, and total 
dissolved inorganic carbon (DIC), w. The wavy line represents the air-sea 
interface, the upper horizontal thick line divides the shallow ocean from the 
deep sea, and the lower horizontal dashed line represents the sediment- 
seawater interface. Concentration fluxes are indicated by unidirectional 
arrows. The sediment-seawater interface is dashed to indicate that there 
is no dynamical distinction between the deep sea and sediments. 
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of CO 2 and CO 2- makes it possible to temporarily forestall a 
discussion of the changing chemical equilibrium. 

The flux ji n represents the rate at which dissolved CaCC >3 is 
added to the oceans by rivers in terms of an equivalent change in 
the whole-ocean concentration of dissolved CaCC> 3 , expressed 
in units of qmol-kg _ 1 -y _1 . The flux B > 0 denotes the rate at 
which CaCC >3 is output from the system, which, in the real ocean, 
corresponds to burial in sediments. However, the model makes 
no dynamical distinction between the deep sea and sediments. 
Thus, B represents the rate at which CaCC >3 exits the shallow 
layer without returning. The dimensionless term v determines 
the strength of an external source of CO 2 , such as volcanic 
emissions. 

The flux R > 0 represents the extent to which upper-ocean 
production of CO 2 by respiration exceeds a theoretical baseline. 
The small change in alkalinity associated with respiration (33) 
is neglected. When v = 0, the baseline case R = 0 is associated 
with a steady state in which w = wo and the riverine influx of dis¬ 
solved CaCC >3 is precisely balanced by burial (i.e., B = jin). In 
this case, the model’s upper ocean would act merely as a way sta¬ 
tion for the temporary storage of dissolved carbonate. Positive 
R, on the other hand, acts to increase w above wo. The model 
assumes that this increase is resisted linearly, so that if the base¬ 
line R = 0 were restored, w would return to wo at the timescale 
t w . At long timescales (3>10 2 y) in the real ocean, the processes 
responsible for such dissipation include the riverine influx of dis¬ 
solved CaCC >3 and the dissolution and precipitation of CaCC >3 
on the seafloor (34, 35). These processes act as the ocean’s 
“homeostat” or “pHstat” at a wide range of timescales (2, 35). 
Aspects of the homeostat that act simultaneously on DIC and 
alkalinity respond to the difference between and B. However, 
the homeostat must also influence DIC independently by oppos¬ 
ing the respiration feedback; otherwise DIC would grow without 
bound when R > 0. Eq. 5 represents this negative feedback by the 
term proportional to w — wo. The homeostat’s dominant charac¬ 
teristic timescale, denoted here by t w , is about 10 ky (36-38) in 
the modern ocean. 

Eqs. 4 and 5 describe the general form of a dynamical sys¬ 
tem in which the riverine influx of alkalinity and DIC produce 
internal feedbacks that collectively act to determine the extent 
to which alkalinity and DIC accumulate before exiting the sys¬ 
tem. The description of the model is completed by the functional 
specification of the feedbacks B and R, to which we now turn. 

The Burial and Respiration Fluxes. The export and burial of CaCC >3 
in the real ocean depends on the physical environment (e.g., 
temperature and pressure), the chemical environment, and eco¬ 
logical factors affecting biogenic carbonate precipitation (e.g., 
phytoplankton community composition). When seawater is 
undersaturated with carbonate ions, CaCC >3 dissolves; when it is 
supersaturated, CaCC >3 can precipitate (33). At low saturations, 
physical and biogenic precipitation—and therefore carbonate 
output—is impeded. The dependence of precipitation on pres¬ 
sure manifests itself as a critical depth beneath which carbonate 
dissolves. Based in part on earlier work (39), I express these 
dependencies by 

B(a, w) = B[c(a, w)] = bj in s(c, c p ), [ 6 ] 

where the carbonate ion concentration c(a,w) derives from the 
equilibrium of the carbonate system (33) and 

S (c,c p ) = ^ J [7] 

is the sigmoidal function shown in SI Appendix, Fig. SI. The 
exponent 7 parameterizes the sharpness of the transition, c p 
is the transitional carbonate ion concentration where s = 1 / 2 , 


and bj in denotes the maximum rate of carbonate burial. The 
apparent requirement that c be computed from a and w would 
appear to be a serious complication, but a change of vari¬ 
ables discussed below solves this problem. The introduction 
of the exponent 7 may also appear unwelcome; however, this 
paper’s conclusions require only that the function s be sigmoidal. 
SI Appendix, Fig. SI suggests that 7~4 corresponds to the 
modern ocean. 

Respiration is also sensitive to environmental conditions. For 
example, warmer temperatures should result in increased res¬ 
piration rates relative to production rates (27, 40-42). Possible 
consequences include an upward shift of respiration and a pos¬ 
itive feedback for shallow-ocean and atmospheric CO 2 levels 
(27, 41, 43, 44). Alternatively, carbonate production by pelagic 
calcifiers may decrease at lower pH (45-48). The “ballast hypoth¬ 
esis” (49-51) suggests that the export of organic carbon out 
of the photic zone is facilitated by its association with bio¬ 
genic carbonate. A decreased supply of ballast would decrease 
export of organic carbon and therefore increase respiration rates 
in the shallow ocean. A positive feedback is again possible 
(27, 52, 53). 

Neither the ballast (54-58) nor the temperature (42, 59-61) 
feedback is free of ambiguity. However, both are reasonable 
hypotheses. This paper addresses the case of the ballast feed¬ 
back, with the expectation that the lessons learned will inform 
consideration of other possible respiration feedbacks, including 
those involving temperature. 

The essential idea is that lower concentrations of CO 2 ^ lead to 
less biogenic carbonate production and therefore less ballast, less 
associated organic carbon export, and more upper-ocean respira¬ 
tion (27, 45-53). These observations suggest that the respiration 
flux R decreases with increasing carbonate ion concentration c. 
For c below a low concentration of carbonate ions, R should 
be near its maximum. Conversely, for c above a high concen¬ 
tration, respiration in the shallow layer is assumed to reach 
its baseline R = 0. The baseline case therefore corresponds to 
the upper-ocean respiration flux that would continue to occur 
despite elevated production of carbonate ballast. 

This reasoning suggests that R decreases with c similar to the 
way B increases with c (SI Appendix, Fig. SI), but with a possibly 
different amplitude and transitional concentration. I therefore 
write 

R(a,w) = R[c(a,w)\=9j in s(c, c x ), [8] 

where s = 1 — s, 9j m is the maximum value of R, and c x is the 
“crossover” concentration that marks the halfway point in the 
transition from high to low R. As for the burial function of 
Eq. 6 , the sharpness of the transition parameterized by 7 is a 
minor detail. SI Appendix, Fig. S2 illustrates how the burial and 
respiration fluxes of Eqs. 6 and 8 modify Fig. 2. 

Transformation from a to c. SI Appendix shows that the rate of 
change of the CO 2- concentration can be approximated by 

c=f(c)(a-w). [9] 

Here the “buffer function” 

/(c) =/o-^ [10] 

approximates the magnitude of the partial derivatives dc/da and 
dc/dw via specific values of /o, /3, and c/ related to the equilib¬ 
rium of the carbonate system (SI Appendix, Fig. S3). Substitution 
of Eqs. 6 and 8 into Eqs. 4 and 5 and insertion of the resulting 
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expressions into Eq. 9 then remove the implicit dependence of c 
on a and w. The new system reads 

C//(C) = jLt[l — bs(c, Cp) — 9s(c, Cx) — v] + w — wo [11] 
w = p[l- bs(c, c p ) + 9s(c, Cx) + u] — w + wo, [12] 

where time is nondimensionalized by ft— t/r w and 

P = jin'i~w [12] 

is a characteristic concentration. The system of Eqs. 11 and 12 
is henceforth called the carbon-cycle model. Its parameters are 
listed in SI Appendix, Table SI. Of these, only p, b, 9, Cx, c v , 
and v are of interest. The remaining 5 parameters are set to cor¬ 
respond to the equilibrium chemistry or the properties of the 
modern ocean. Although they are retained in the analysis that 
follows, they serve only to maintain realism. 


Steady State and Its Stability. Steady states, or fixed points, occur 
where c = w = 0. Under this condition, the addition of Eqs. 
11 and 12 yields bs(c, c p ) = 1. Substituting Eq. 7 and solving 
for c, one then finds the unique steady-state carbonate ion 
concentration 

c* = (b— l) - 1 ^' y Cp, [14] 


which exists only for b > 1. Substitution of c* for c in Eq. 12 then 
provides the steady-state DIC concentration 

w* =w 0 +p\9+ v - —— gCp 7 — 7 ). [15] 

\ (b 1) Cx T Cp J 

Note that the untransformed model (Eqs. 4 and 5) yields the 
same steady state after substitution of Eqs. 6 and 8 . SI Appendix 
shows that parameter values consistent with the modern ocean 
(SI Appendix, Table SI) predict values of c* and w* that are also 
consistent with the modern ocean. 

SI Appendix also derives and interprets the results of a linear 
stability analysis. The requirements for stability are expressed 
in terms of the trace Tr and determinant A of the Jacobian 
matrix of Eqs. 11 and 12 evaluated at the fixed point. These 
calculations yield 


Tr = —1 — 


A 

2 ~ 



9bc^c2 \ 

[(&-i)cg + c?]V’ 


[16] 


where 


A = 


2M7/o(&-1) 




[17] 


b(b - 1)0/7 cjf + 6 c£ 

The fixed point is unstable when Tr > 0. Fig. 3 provides sta¬ 
bility boundaries in the planes of p and b, p and c x , and p 


and 9. Instability is favored for large p=ji R r w because greater 
rates of riverine input and longer damping times provide for 
faster-growing fluctuations. For some combinations of parame¬ 
ters, however, homeostatic damping always maintains stability. 

Limit Cycle 

SI Appendix shows that the fixed point becomes unstable via a 
Hopf bifurcation, which leads to periodic nonlinear oscillations 
called limit cycles. Fig. 4 depicts a stable limit cycle in the phase 
plane of c and w. SI Appendix, Fig. S4 illustrates the time series 
c(t) and w(t) on the limit cycle, along with the time evolution of 
quantities derived from these 2 quantities. 

To understand how the limit cycle operates, consider the evo¬ 
lution of the periodic trajectory from the point where w is at 
its minimum, i.e., at the lower crossing of the nullcline w = 0 . 
Here, the ballast feedback begins to increase DIC by the addi¬ 
tion of CO 2 , which acts to consume CO3 “ ions via the reaction 
of SI Appendix, Eq. S 8 . However, the relative rate at which 
the CO 3 concentration decreases compared with the rate at 
which DIC increases becomes progressively slower as the ocean 
becomes less buffered at lower CO 3 concentrations. Eventually 
the CO 3 - concentration reaches a minimum where the loss of 
CC> 3 _ via burial and the CO 2 feedback is balanced by the home¬ 
ostatic response, and it crosses the nullcline c = 0. However, DIC 
continues to rise until it is overcome by the homeostatic feed¬ 
back, reaching a maximum where it crosses the nullcline w = 0 . 
The homeostatic feedback then acts to decrease w while increas¬ 
ing c. The CC> 3 _ concentration reaches its maximum when 
carbonate burial becomes stronger than the addition of dissolved 
CaCC >3 from rivers and the weakened homeostatic feedback. A 
rapid upward jump of the burial rate then acts to reduce both 
the CO 3 and DIC concentrations. Eventually DIC reaches its 
minimum and the cycle repeats. 

A single period of the limit cycle embodies several attributes 
of the marine carbonate system that have attracted consider¬ 
able attention in recent years. The rise of DIC, for example, is 
initially accompanied by a rise of CO 2 (SI Appendix, Fig. S4C) 
and a reduction of pH (SI Appendix, Fig. S4 E). This process is 
called ocean acidification (62). The acid is eventually neutralized 
by the dissolution of seafloor carbonate and other homeostatic 
processes, but only slowly, at a modern timescale of about 10 
ky—part of the long legacy of any rise in CO 2 levels (34). Finally, 
the surfeit of dissolved CaCC >3 must eventually be buried. In 
the limit cycle pictured here, carbonate burial occurs as a sharp 
pulse (SI Appendix, Fig. S4 H) analogous to the “cap-carbonate” 
deposits that follow global glaciations (63). 

Periodic disruptions of the carbon cycle independent of peri¬ 
odic forcing may have occurred in the geologic past (64-67). 
The following section shows, however, that the potential exis¬ 
tence of a limit cycle is more important than its actual existence. 
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Fig. 3. (/4-C) Stability in the planes of p and b (A), and c x (6), and /x and 6 (C). Values of fixed parameters are given in SI Appendix, Table SI. Crossing the 
red and blue boundaries results in supercritical and subcritical Hopf bifurcations, respectively. The fixed point is stable only where indicated. However, the 
region of stability of the limit cycle extends into the region of the stable fixed point adjacent to the subcritical bifurcation. Fig. 6 maps the bistable region 
in the plane of c x and 0. 
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Fig. 4. A stable limit cycle in the phase plane of the CO^ concentration c 
and the DIC concentration w. Values of parameters are given in SI Appendix, 
Table SI. The dotted curves represent the nullclines c = 0 and w = 0. 


Three aspects of the limit cycle are central to that discussion: Its 
amplitude, its period, and its dynamical origin. 

Unlike linear oscillations, the amplitude or size of a limit 
cycle is independent of the initial condition. As shown in Fig. 4, 
trajectories originating inside the limit cycle grow outward, 
whereas trajectories originating outside move inward. The result¬ 
ing amplitude is therefore a property of the system itself rather 
than its perturbation. 

However, the size of the limit cycle does depend on the sys¬ 
tem’s parameters. An appropriate measure of size along the DIC 
axis is the difference Aw — w max — w*, where w max is the maxi¬ 
mum value of w in the cycle; one similarly has Ac = c max — c*. 
SI Appendix shows that Aw and Ac generally increase with p 
and 6 but decrease with b. The limit cycle’s size grows roughly 
linearly with p because p sets the rate of all processes other 
than homeostatic damping. Because rates and size scale similarly 
with p, the dimensionless period T of all but the smallest limit 
cycles should be effectively independent of p. These observations 
imply that 

T~T 0 + 2A/p, [18] 

where A = (AwAc) 1 / 2 is a measure of the amplitude of the 
limit cycle in units of pmol-kg -1 and To ~ 1 is the dimension¬ 
less period of an infinitesimal limit cycle. The factor of 2 reflects 
an assumption that roughly twice as much time is required to 
go through a complete cycle compared with the time it takes to 
change w or c by its typical excursion. 

SI Appendix, Fig. S5 shows that Eq. 18 is a reasonable approx¬ 
imation. It characterizes the limit cycle well partly because the 
complicated dependence of A and T on the system’s parameters 
other than p (SIAppendix) is similarly expressed by A and T. The 
physical interpretation of Eq. 18 is readily exposed by rewriting 
it in terms of dimensional time. Setting r = Tt w / 2, to = Tqt w /2, 
and recalling p=j m r w , one finds 

A — in(T — to). [19] 

This relationship, which essentially follows from dimensional 
analysis, shows that, for t^>t o, the limit cycle’s size is roughly 
given by the riverine influx jj n integrated over the time required 
for the system to be forced maximally out of its steady state. 

The manner in which the limit cycle arises is also of interest. 
There are 2 kinds of Hopf bifurcations—subcritical and super¬ 
critical (21, 28, 29). SI Appendix performs a standard calculation 


(21, 28, 29) to determine where each Hopf bifurcation occurs 
here; the results are color coded along the curves of Fig. 3. In 
the carbon-cycle model, the subcritical Hopf bifurcation is adja¬ 
cent to a bistable region where the stable fixed point and a stable 
limit cycle appear with an unstable limit cycle. The following 
section establishes the importance of this region for the inter¬ 
pretation of the model’s response to perturbations of the stable 
fixed point. 

Excitability 

Now suppose the system is in a stable steady state with no 
injection of CO 2 (v = 0). At time t = 0 we set iz = 0.35, which 
corresponds to injecting additional DIC into the oceans at a rate 
equal to 35% of riverine inputs. From Eqs. 14 and 15, we see 
that CO 2 injection increases w * by i/p, but c* is unchanged. 
The new fixed point retains its stability, because stability is 
independent of v (Eqs. 16 and 17). However, the state of the 
system must now relax to the new fixed point. Fig. 5 A and B 
shows that there is a small excursion beyond the fixed point 
as the trajectory spirals inward to the new steady-state DIC 
concentration. 

Next we repeat the same procedure, but with a slightly larger 
CO 2 injection rate, i/ = 0A0. The initial condition remains the 
same and the new, stable steady-state DIC concentration is only 
marginally higher than before. However, Fig. 5 C and D shows 
that the excursion beyond the fixed point is now roughly an order 
of magnitude greater than before. Indeed, its size is compara¬ 
ble to the large limit cycle of Fig. 4. The trajectory, however, 
is not periodic, and the system soon relaxes back to its stable 
steady state. 

SI Appendix, Fig. S 6 shows how the size of such excursions 
varies with u. The size is defined by the difference between 
the maximum values of w and w*, where w* is the DIC fixed 
point computed from Eq. 15 for a given value of v. SI Appendix, 
Fig. S 6 shows that the transition from small to large excursions is 
sharp yet continuous. Here it occurs near v = 0.391. Beyond this 
threshold, the size of the excursion no longer grows appreciably 
with v. 

Dynamical systems that exhibit phenomena similar to those of 
Fig. 5 and SI Appendix, Fig. S 6 are called excitable (21). When a 
stimulus such as v exceeds a threshold, excitable systems produce 
a transient excitation significantly larger than the stimulus, as 
seen here. Excitable systems were first introduced to understand 
the behavior of action potentials (nerve impulses or spiking) in 
neurons ( 68 , 69). They have since been studied in the context of 
numerous other problems, including glacial cycles (32) and the 
response of peatlands to warming (70). 

Excitations occur near Hopf bifurcations (21). Even though 
the fixed point is stable, a remnant of the stable limit cycle 
exists in the phase plane. Here, as illustrated in Fig. 5, when v 
is above a critical threshold v c , the trajectory effectively hops 
on to a vestige of the limit cycle before settling in to the 
fixed point. The size and timescale of excitations are there¬ 
fore comparable to limit cycles. Specifically, their amplitudes 
are insensitive to the extent to which the threshold parameter 
exceeds the threshold (SI Appendix, Fig. S 6 ), and their dura¬ 
tion is comparable to the period of limit cycles. Unlike limit 
cycles, however, they do not require that a parameter of the 
system, such as p, be increased beyond a point of bifurcation. 
For the carbon-cycle model studied here, this means that excita¬ 
tions can occur in a stable carbon cycle that is merely subjected 
to an above-threshold injection of CO 2 while its constituent 
parameters—e.g., riverine fluxes into the system, burial fluxes 
out of it, and the timescale of homeostatic equilibration—remain 
the same. 

Fig. 6 shows where excitations occur and at what value of v c , 
in the plane of c x and 6, along with the stability boundary at 
which the Hopf bifurcation occurs. Fig. 6 also shows the bistable 
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Fig. 5. Evolution following perturbations below and above the excitations threshold. ( A-D ) Phase-plane representation and the time series w(t) for a sub¬ 
threshold perturbation with v = 0.35 (A and 6) and an above-threshold perturbation with v = 0.40 (C and D). The initial condition and all other parameters 
are the same in each simulation. The crossover CO^' concentration c x = 55 pmol-kg -1 ; other parameters are given in SI Appendix, Table SI. The larger 
value of v in the above-threshold case increases the DIC fixed point w* by 0.05 /j = 12.5 pmol kg -1 in C compared with A, which is nearly imperceptible in 
the plots. 


region adjacent to the stability boundary where the subcritical 
Hopf bifurcation occurs. One sees that excitations occur in a wide 
region of parameter space adjacent to the bistable region. 

The relation between excitations and bistability is further clar¬ 
ified in Fig. 7, which provides a ID view along the c x axis. The 
bistable region occurs at values of c x where the stable fixed point 
and stable limit cycle coexist with an unstable limit cycle. In this 
region, values of v that depress the minimum w of the stable limit 
cycle below the minimum w of the unstable limit cycle create 
an immediate jump to the stable limit cycle. At the left bound¬ 
ary of the bistable region, the unstable limit cycle and the stable 
limit cycle collide via a saddle-node bifurcation of cycles (29). 
For values of Cx below this point, only the fixed point is stable. 
Whereas an above-threshold v creates a jump to the stable limit 
cycle when c x is in the bistable region, a similar above-threshold 
v creates an excitation at smaller values of c x . Fig. 8 illustrates 
how the unstable limit cycle determines the evolution of phase- 
plane trajectories in the bistable regime. SI Appendix, Fig. S7 
shows that excitations near the saddle-node bifurcation of cycles 
can create a finite train of pulses before returning to the stable 
fixed point. 

Discussion 

If the carbon cycle is indeed excitable, the great events of the 
geologic past should reveal signs of such dynamics. Among these 
events, those associated with mass extinction are of particu¬ 
lar interest. It is also of interest to ask how the hypothesis of 
excitability informs our understanding of the risks of human- 
induced carbon-cycle disruptions in the future. The following 
discussion provides some first steps in these directions. 

Characteristic Events. Because excitations inherit the properties 
of limit cycles, Eq. 19 relates not only the amplitude and 
period of limit cycles via the flux j\ n , but also the amplitude 
and duration of excitations. As expected, numerical simula¬ 
tions of the carbon-cycle model show that the excitation size 
Aw)~7i n T for large r, where r is the time from the onset of 
an excitation to its peak {SI Appendix, Fig. S8). The geochemi¬ 
cal record provides data from which the normalized size A w/w* 
may be estimated for real events; it also provides estimates of 
r (7). If real events are excitations, then the above reason¬ 


ing suggests that their size and timescale should approximately 
satisfy 


where a is a constant. Because aj m /w* is a concentration flux 
divided by a steady-state concentration, I call it the specific rate 
of the excitation. 

Fig. 9 plots Aw/w* as a function of r for 31 global disrup¬ 
tions of the carbon cycle during the last 542 My (7). Five of the 
events are associated with major mass extinctions. Variations in 
t are likely related to changes in the timescale of homeostatic 
damping over geologic time (7). The straight line and its uncer¬ 
tainty (the yellow region) predict the size of the disruption that 
would occur if it were quantitatively related to the cessation of 
organic carbon burial and the resulting accumulation of CO 2 (7). 
Roughly half of the events occur in this region; they each display 
surges of carbon that grow with approximately the same specific 
rate. In addition to being well populated, the yellow region also 



25 50 75 100 125 150 

Cx [/rmol kg' 1 ] 


Fig. 6. Contours of v c (dotted) in the region of excitations (pale yellow) 
and jumps to the stable limit cycle in the region of bistability (pale blue) in 
the plane of c x and 0. The blue and red portions of the stability boundary 
correspond to a subcritical and a supercritical Hopf bifurcation, respectively. 
Fixed parameters are specified in SI Appendix, Table SI. Excitations or jumps 
are defined to occur if there exists a v for which dA w/dv> 10 3 , where 
Aw = w max — w*. The threshold v c is the value of v where that derivative is 
greatest. 
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Fig. 7. Bifurcation diagram for the parameter c x , the crossover CO^ 
concentration for the ballast feedback. The solid and dotted black line 
respectively represents the stable and unstable fixed point w*. The solid 
blue line indicates the maximum and minimum values of w in the stable limit 
cycle and the dashed red line represents the same extremes of the unstable 
limit cycle. The subcritical Hopf bifurcation occurs where the radius of the 
unstable limit cycle goes to zero, at c x = 62.61 pmol-kg -1 . The saddle-node 
bifurcation of cycles occurs where the unstable and stable limit cycles col¬ 
lide, at c, = 55.89 u mol - kg 1 Excitations occur at smaller values of c x , to 
the left of the arrowhead. Fixed parameters are specified in SI Appendix, 
Table SI. 


thing in common other than their slowness. They likely rep¬ 
resent a quasistatic (10) response to subthreshold changes in 
the parameters of the carbon cycle, such as changes in burial 
rates. These minor events therefore probably embody the clas¬ 
sic proportionate response. The anomalously fast mass extinc¬ 
tion events in the upper region should, however, be related to 
excitability. To understand how, I first discuss the excitation 
threshold. 

The Excitation Threshold. Practically speaking, the amplitude of 
an excitation should significantly exceed the amplitude of the 
perturbation required to exceed the threshold (21). In terms of 
the carbon-cycle model, such a definition immediately provides 
a conservative upper bound for v c by requiring that the increase 
in w* due to v be smaller than the excitation amplitude Aw. 
Obtaining the former from Eq. 15 and the latter from Eq. 20 , 
one then has v c p, < aj m r with a ~ 1. Similar reasoning applies 
to the real carbon cycle. Given a damping timescale t w , the DIC 
fixed point should again increase by roughly u c p, at the threshold. 
The assumption that characteristic events are excitations then 
provides the same bound on this increase, but with a~0.1. 
In either case, the excitation timescale r should be similar 
to the damping timescale. Recalling pi = ji n r w , one then finds 
v c < q or 


appears to separate benign events from 4 of the 5 mass extinc¬ 
tion events. The disruptions within the yellow region therefore 
appear to share important characteristics in addition to their 
common specific rate. Henceforth these disruptions are called 
characteristic events. 

I hypothesize that characteristic events represent excitations 
of the marine carbon cycle. If so, then the size of characteristic 
events is predicted by Eq. 20 , just as it is for model excita¬ 
tions. The yellow region in Fig. 9 satisfies Eq. 20 for a ~ 0.1 (7). 
Model excitations scale similarly, but with a ~ 1 (SI Appendix, 
Fig. S8). This difference likely reflects mechanisms, such as 
limitations imposed by organic carbon burial (7), that are not 
captured by the model. The important point, however, is that for 
both observed disruptions and model excitations, the amplitude 

Aw OCJjnT. 

The shared scaling of characteristic events and model excita¬ 
tions with the riverine flux j m points to a common foundation. 
The carbon-cycle model is an open, nonlinear chemical reac¬ 
tor driven out of equilibrium by a flux of reactants into and out 
of the system. So, too, is the real ocean; elemental concentra¬ 
tions are determined by the balance of inputs and outputs while 
buffering is controlled by heterogeneous carbonate equilibria 

(71) . More generally, instability, bistability, and oscillations can 
occur in nonlinear chemical systems forced out of equilibrium 

(72) . In such problems, the size of limit cycles—and therefore 
the size of any excitations—is a property of the system rather 
than its perturbation (29). The model expressed by Eqs. 11 and 
12 shows how this works in an idealized carbon cycle. The pro¬ 
portionality of event sizes to the driving flux is a robust property 
of the model that should be equally applicable in more com¬ 
plex settings. The common scaling with j m also provides strong 
evidence that characteristic events display an intrinsic nonlinear 
response to perturbation rather than a proportionate response 
to an external stress. Because different above-threshold stressors 
would yield a similar response in an excitable system, the exci¬ 
tation hypothesis also explains why the yellow region in Fig. 9 is 
relatively well populated. In addition, because the carbon iso¬ 
topic signals typically return to their initial steady state after 
each event (e.g., Fig. 1), the steady state apparently remains sta¬ 
ble. The combination of these observations is consistent with 
excitability. 

Events lying below the yellow region in Fig. 9 occur at 
relatively slow specific rates and do not appear to share any¬ 


Vc 



model 
real world. 


[ 21 ] 

[ 22 ] 


Fig. 6 shows that the excitation threshold in the carbon-cycle 
model is consistent with Eq. 21 near the subcritical Hopf 
bifurcation. 

The combined output of CO 2 from modern subaerial and sub¬ 
marine volcanism produces about 0.06 Pg-C-y -1 (73), which is 
about 15% of the modern riverine flux of 0.41 Pg C-y -1 (74). If 
the volcanic flux were doubled, the perturbation expressed by v 
would therefore be about 0.15, equivalent to the predicted upper 
bound for v c (Eq. 22) and therefore enough, in principle, to 
initiate an excitation. 


Mass Extinctions. Fig. 9 shows that 4 mass extinction events 
exhibit a specific rate significantly greater than that of 
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Fig. 8. Phase-space trajectories in the bistable regime. Here c x = 
58 (rmol kg -1 and all other parameters are the same as in Fig. 5 C and 
D. The red dashed limit cycle is unstable; thus trajectories initialized inside it 
return to the (stable) fixed point. Trajectories initialized outside the unstable 
limit cycle evolve to the stable limit cycle. 
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Fig. 9. The relationship between the relative size and duration of 31 
disruptions of the global carbon cycle during the last 542 My (7). The rel¬ 
ative size Aw/w* is obtained from changes in the isotopic composition 
of carbonate carbon that occur in time series similar to that of Fig. 1. 
The yellow region contains characteristic events; these events satisfy Eq. 
20 with a ~ 0.1. The events in the pale red region, which include 4 of the 
5 mass extinction events, grow faster than characteristic events. The light 
blue region contains minor events with relatively slow growth rates. The 
labeled events are associated with the end-Cretaceous (KT), end-Triassic (TJ), 
end-Permian (PT), end-Ordovician (Ord), and Frasnian-Famennian (FF) mass 
extinctions. Ref. 7 provides descriptions of each event and a discussion of 
uncertainties. 

characteristic events, by a factor of 2 to 5 (7). These events 
can be characterized as excitations driven by an injection rate 
v significantly greater than v c . 

The idea follows reasoning similar to that used to bound v c . 
Excitations at threshold exhibit the characteristic size Aw ~ 
aj m r (Eq. 20 ). Suppose that CO 2 injection occurs at rate ( v c + 
e)jin, where e > 0. Eq. 15 then shows that Aw will increase by efi 
relative to its size at the threshold, resulting in 

Aw ~in(ar + £T W ) . [ 23 ] 

In other words, above-threshold excitations are larger than 
threshold excitations by a factor of about 1 + ST w /ar. Since 
the excitation timescale r should not appreciably change, the 
observed rate would similarly increase. Because t ~ t w , the main 
interest lies in e/a, which is on the order of 10 e for the real car¬ 
bon cycle. So, for v c ~ 0.1 (the upper bound of Eq. 22 ), excitation 
at twice the threshold would increase Aw and its growth rate by 
roughly a factor of 2 . 

The association of massive flood-basalt volcanism with the 
end-Permian (23), end-Triassic (22), and end-Cretaceous (24, 
25) mass extinctions suggests how this would work. Flood basalts 
release CO 2 at a rate of about 3.5 x 10 -3 Pg C/ km 3 of magma 
(75). Schoene et a!. (24) estimate that up to 11 km 3 -y _1 of 
magma erupted from Deccan volcanism over a 10-ky period 
tens of thousands of years before the end-Cretaceous extinc¬ 
tion. The resulting CO 2 flux is equivalent to about 10% of the 
modern riverine flux, yielding ^~ 0 . 1 , commensurate with the 
upper bound for v c . Other episodes of massive volcanism (22- 
25) may exhibit similar pulses (76). Such events may indeed be 
well above threshold: If v c were much smaller than its upper 
bound, then e = 1 / — v c ~ 0.1, yielding a/e ~ 1. The size and rate 
of the event would then be roughly double those of characteristic 
disruptions. 

Excitation of the Modern Carbon Cycle. As a result of anthro¬ 
pogenic CO 2 emissions, Earth’s oceans currently absorb about 
2.3 Pg C-y -1 (77), which is about 5.6 times the riverine carbon 
influx. The resulting excitation parameter, v = 5.6, is more than 
50 times greater than the predicted upper bound for v c . That 
prediction, however, assumes a continuous, constant injection of 
CO 2 . If CO 2 injection occurs instead for a time much less than 


the damping timescale, its ability to initiate an excitation may be 
attenuated by the ocean’s homeostatic response. 

Whether or not excitation ensues will nevertheless still depend 
on the injection rate. Suppose that CO 2 injection occurs over a 
time U <t w . If excitations follow when v>v c for U>t w , then 
excitations should still occur when U < t w if the DIC addition 
uji n ti is approximately equal to or greater than the total input 
resulting from injection at the threshold v c over 1 damping 
time. In other words, there is a t^-dependent threshold v' c that 
approximately satisfies v' c ti = v c t w for small U. Therefore, 

v c (t i )~i/ c T A, U<t w . [ 24 ] 

ti 

The scaling like t” 1 should increase in accuracy as ti—>-0. SI 
Appendix, Fig. S9 shows that numerical solutions of the carbon- 
cycle model compare well to this prediction. The dependence of 
the short-time threshold on U is qualitatively similar to the notion 
of a critical ramping rate (70). 

Given the real-ocean upper bound v c <0.1 (Eq. 22 ), sub¬ 
stitution of t w = 10 4 y and ti = 10 2 y into Eq. 23 yields the 
century-scale upper bound v' c < 10. Mean 21st-century oceanic 
carbon uptake will likely range from 1.8 to 4.2 PgCy -1 
(77), which translates to 4.4 < v < 10.2, similar to the upper 
bound for v' c . With respect to the excitation threshold, the 
modern perturbation is therefore approximately equivalent 
to the aforementioned end-Cretaceous pulse; its duration is 
2 orders of magnitude shorter but its flux is 2 orders of 
magnitude greater. Fig. 10 illustrates this correspondence 
graphically. 

The short-time critical rate v' c may be alternatively expressed 
as a critical mass m c . Multiplication of both sides of Eq. 24 by 
jinti shows that injection with v>v' c over a timescale ti < t w is 
equivalent to adding a mass of carbon greater than 

m c ~ 12u c j m r w V, [ 25 ] 

where V = 1.35 x 10 21 kg is the mass of the oceans (2) and 
the factor of 12 converts mol to g. Ref. 7 derives an equiva¬ 
lent expression from different assumptions. In conjunction with 



Fig. 10. Equivalence, with respect to the excitation threshold, of the per¬ 
turbations of the modern and end-Cretaceous carbon cycles, expressed in 
terms of the dimensionless C0 2 injection rate v. The straight line labeled 
1 /' is the threshold's upper bound predicted by Eq. 24, assuming the mod¬ 
ern damping timescale r w = 10 4 y; the segment labeled z/ c (Eq. 22) provides 
the upper bound for times t; > 10 4 y. The circles represent projected 21st- 
century ocean carbon uptake rates for a range of plausible scenarios (77); 
the vertical line indicates their uncertainty. The symbols labeled KT pro¬ 
vide the median and upper limit of the 68% confidence interval estimated 
for the period of peak Deccan volcanism tens of thousands of years before 
the end-Cretaceous extinction (24). Because the excitation threshold scales 
like 1 /ti, both the modern and end-Cretaceous perturbations potentially lie 
near the threshold's upper bound despite the 2 orders of magnitude that 
separate their rates. 
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Eq. 22 , Eq. 25 predicts the upper bound m c < 405 Pg C, which 
falls within the projected uptake of 290 to 540 Pg C from the year 
1850 to 2100 (77). Thus, in terms of either 21st-century uptake 
rates or total postindustrial uptake, the modern marine carbon 
cycle may be at the precipice of excitation. 

Finally, there is an important caveat. Since the carbon-cycle 
model applies only to timescales greater than that of ocean 
mixing [about 1,000 y (2, 3)], it does not necessarily inform 
understanding of phenomena at timescales shorter than 1,000 y. 
However, if the addition of CO 2 to the oceans were not apprecia¬ 
bly damped at submillenial timescales by processes not discussed 
here, the reasoning outlined above would remain valid at these 
short timescales. 

Conclusion 

Major disruptions of Earth’s carbon cycle are typically inter¬ 
preted as a proportionate response to an environmental per¬ 
turbation ( 8 , 16, 20). This paper suggests instead that many 
of these events represent the nonlinear amplification of pro¬ 
cesses that operate within the carbon cycle. To illustrate how 
this could work, this paper constructs and studies a simple 
dynamical-system model of the marine carbon cycle. The results 
suggest that disruptions can follow when a stable steady state 
is perturbed beyond a threshold. A likely source of pertur¬ 
bations is an increased flux of CO 2 into the oceans. The 
autocatalytic amplification of model disruptions exhibits a char¬ 
acteristic increase and rate of growth of the ocean’s store 
of carbon before returning to the original steady state. Such 
characteristic sizes and rates, independent of the history of 
external perturbations, are classical features of nonlinear sys¬ 
tems (29). The historical record of the real carbon cycle also 
exhibits similar characteristics (7). These shared properties sug¬ 
gest that the model displays important features of the real 
carbon cycle. 

These findings suggest that a characteristic excitation of upper- 
ocean CO 2 levels follows the injection of CO 2 into the oceans 
at an above-threshold rate. If injection continues over a time 
greater than the timescale t w over which the oceans homeostat- 
ically adjust to changes in pH, the threshold’s upper bound is 
roughly equivalent to the higher estimates of the rate at which 
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CO 2 degasses from major flood-basalt eruptions associated with 
mass extinction (24, 76). If injection is instead limited to a shorter 
duration U, the threshold increases by a factor of t w /U. Thus, the 
relatively slow rate of CO 2 injection commensurate with massive 
volcanism at geologic timescales turns out to be roughly equiva¬ 
lent, in terms of its potential to reach the threshold, to the much 
stronger but briefer perturbation of the modern carbon cycle. 
Mass extinction events appear to be associated with excitations 
well above threshold. 

These conclusions rely in part on the assumption that the 
dynamics of a 2D dynamical system represent a subset of 
the complex behavior exhibited by the real carbon cycle. This 
assumption does not require that the mechanisms encoded by 
the terms within the dynamical system be strictly correct. How¬ 
ever, it does demand that the model’s qualitative dynamical 
properties, such as the possible existence of a limit cycle, be 
realistic. This may not be true. One could argue, for example, 
that phenomena at submillenial timescales, such as ecological 
reorganization, will act as strong negative feedbacks in local envi¬ 
ronments, impeding global autocatalytic self-organization. The 
present state of the geochemical record makes such notions dif¬ 
ficult to rule out observationally. On the other hand, numerical 
simulation of more detailed carbon-cycle models (44) should 
allow one to test whether excitations are expressed in such frame¬ 
works. Likewise, careful analysis of individual disruptions or 
specific periods in the geochemical record should allow one to 
test whether quantitative signatures of excitable systems exist. 
The possibility of resonance with the fluctuations of external 
perturbations (21), such as those induced by orbital variations 
(12), offers one such route. The work reported here identifies 
why the carbon cycle may be excitable, how excitability may 
have been expressed in the past, and why an excitation may 
occur in the future. In a curious twist of geological and biolog¬ 
ical evolution, dynamical mechanisms underlying environmental 
upheaval and mass extinction may be similar to those that make 
neurons spike. 
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Supporting Information Text 

Transformation from a to c via the buffer function /(c). Substitution of Eq. (6) and Eq. (8) into Eq. (4) and Eq. (5) yields 


a = 2j in [l — bs(c, c p )} [SI] 

w = jin [1 - bs(c,c p ) + 0s(c, Cv) +v]— {w - Wo)/t w . [S2] 

Figure S2 illustrates this system schematically. 

The dependence of c on a and w occurs through complex nonlinear expressions of the carbonate system’s thermodynamic equilibrium (1). 
However, c can be obtained from a and w from 

dc _ da dc dw dc 
df df da + dt dw 

Evaluating the partial derivatives—known in different form as buffer factors (2)— also requires knowledge of c(a, w). However a reasonable 
approximation exists. In general, the partial derivatives have opposite signs (2): increasing alkalinity for constant DIC shifts the equilibrium 
toward a greater concentration of carbonate ions, whereas increasing DIC for constant alkalinity moves the equilibrium toward a lower 
carbonate-ion concentration. Figure S3 shows that, as a function of c for a fixed value of w, the magnitudes of the partial derivatives are nearly 
equal. Moreover they can be approximated by another sigmoidal function. In summary, 


dc ^ dc ^ cP 

da ~ dw ~ 0 C P + Cj 


/(c)- 


[S4] 


Numerical solutions of the carbonate equilibrium show that the positive prefactor /o, the exponent /3, and the transition concentration c/ vary 
only weakly with w. Here these parameters are taken to be constant and the approximations in Eq. (S4) are treated as equalities. Then 


c = /(c) ( a-w ). 


[S5] 


Substitution of Eq. (SI) and Eq. (S2) then provides c as a function of c and w. The resulting dynamical system for c and w, identical to Eq. (11) 
and Eq. (12) of the main text, is 


c//(c) = ffl — bs(c,c p ) — 9s(c,c x ) — v\ + w — wo [S 6 ] 

w = ffl — bs(c, c p ) + 9s(c, c x ) + u] — w + wo, [S7] 

where fi = ji n r w and time is expressed in units of r w . 

The buffer function /(c) parameterizes the equilibrium chemistry of carbonate system. While the fit of Figure S3 provides parameters 
corresponding to the true equilibrium, from a dynamical point of view any sigmoidal shape would yield qualitatively similar results. To 
understand this from a chemical perspective, add the first and second reactions in Eq. (3) to obtain 

C0 2 + CO 3 - +H 2 O ^ 2 HCO 3 . [S 8 ] 

If, as in the modem ocean, the concentration of COg - is much greater than the concentration of CO2, the addition of one mole of CO2 (via, 
e.g., the term v) consumes nearly as much CO3 - to form additional HCO^. In this case, the consumption of CO3 - is effectively limited 
by the concentration of CO2 and /(c) is near its maximum value. If we imagine instead that the equilibrium COg - concentration is much 
lower than that of CO2 for the same concentration of DIC, the consumption of CO3 - is effectively limited by the supply of COg - and /(c) is 
much smaller, approaching zero as c —» 0. The sigmoidal shape of /(c) arises from these differences at the limits of small and large c, which 
correspond, respectively, to poorly buffered and well-buffered seawater. 

Characterization of the Steady State. As one would expect, the steady state CO 3 - concentration c* (Eq. (14)) is determined solely by 
the balance between inputs and outputs, and therefore only the parameters b, c p , and 7 that characterize the burial flux given by Eq. ( 6 ). 
Figure S10A provides contours of c* as a function of of c p and b. The average modem deep COg - concentration is 86 /rmol kg -1 (3). Taking 
that as the value for c* and setting 7 = 4 and c p = 110 /rmol kg -1 as in Figure SI, one finds b = 3.7. The parameter b is the ratio of the 
maximum rate of carbonate burial to the rate of input of dissolved CaCOa. It should be similar to the carbonate “overproduction ratio” (4). The 
overproduction ratio indicates how many times a mole of CaCC >3 cycles between biogenic precipitation and dissolution. Broecker (5) estimates 
its value to be 4, which is remarkably consistent with the value of 3.7. 

The fixed point w * (Eq. (15)) depends on several additional parameters. Among these, the characteristic concentration p = ji n T w may 
be estimated for the modern ocean. The modern riverine dissolved influx / n is about 0.025 /rmol kg -1 yr -1 ( 6 ). As discussed in the main 
text, identification of r w with the modern homeostat timescale implies t w ~ 10 Kyr (7-9). Therefore the modern value of fi is about 250 
/rmol kg -1 . The characteristic DIC concentration Wq represents surface DIC in the absence of a ballast feedback or CO 2 injection; the modem 
surface DIC concentration of roughly 2000 /rmol kg -1 (3) is a reasonable choice. The parameters 9 and c x are not straightforwardly measured; 
they instead represent the combined effects of many biogeochemical processes. Plausible values may be inferred, however, from Figure S10B, 
a plot of w* as a function of c x and 9. The modern deep DIC concentration is about 2280 /rmol kg -1 , with about 2100 /rmol kg -1 at a depth 
of about 100 m. Any combinations of c x and 9 that fall between the contours for 2100 and 2280 /rmol kg -1 in Figure S10B are therefore 
consistent with the modem ocean. For example, if c x = 55 /.tmol kg -1 , then 9 would lie between 2.5 and 7.1 
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Linear stability analysis. This section evaluates the stability properties of the fixed point c* ,w* given by Eq. (14) and Eq. (15) of the main 
text. It begins by expressing the calculation in its standard form (10-12). 

First, shift the fixed point to the origin via the change of variables 


v = c — c and u = w — w . 

Then rewrite Eqs. (11-12) in terms of v and u to obtain the system 

v = F(v,u) 
u = G(v,u). 


Next, form the Jacobian matrix 


where the partial derivatives 


J = 


F v F u 
G v G u 


_<TF _ cAF _ dG _dG 

P V 0 1 F U rs > 'GTV rv 5 It rv 

OV OU OV OU 

are evaluated at the origin (v = 0, u = 0). Then compute the trace 

t = tr J = F v + G„ 


[S9] 

[510] 

[511] 

[512] 

[513] 

[514] 


and the determinant 

A = det J = F V G U - F U G V . [S15] 


The stability of the fixed point depends on the sign of the real part of the eigenvalues A of the matrix J — XI, where I is the identity matrix. 
The eigenvalues are given by 


A = 


r ± \/t 2 — 4 A 


It suffices to determine the sign of r: when r < 0, the fixed point is stable, and when r > 0 it is unstable. 
Application of the above recipe yields 


r = -l- ^ [ 1 — 


Obc^cl 


[(6- l)cZ + cl} 2 


where the determinant 


A= 2 M7 /o(&-1)T+ 1 c£- 1 


[S16] 


[S17] 


[SI 8] 


b(b — VfPhcTf + bcp 

Because b > 1 and all other quantities are positive, A > 0, and so too is the term proportional to 6 in Eq. (S17). One then sees that the 
condition for instability, r > 0, can be met, for example, by choosing 6 sufficiently large so that the term in parentheses in Eq. (S17) is 
negative, and choosing fi sufficiently large so that r > 0. Define /x c to be the value of /r where r = 0. At /r = fi c , the eigenvalues (Eq. S 16) 
are pure imaginary: 

A = ±«A 1/2 , [S19] 


where i 2 = — 1. Increasing fi > /j, c then makes real part of A positive and the fixed point becomes unstable. The same conclusion, via a similar 
argument, follows if 6 is treated as the bifurcation parameter, with 6 C its value at the point of bifurcation. 

When the eigenvalues of a two-dimensional dynamical system cross the imaginary axis as they do here, the system undergoes a Hopf 
bifurcation. Determination of whether the Hopf bifurcation is supercritical or subcritical follows the procedure given, for example, on pages 
169-170 of Izhikevich (12), page 289 of Strogatz (11), or pages 152-156 of Guckenheimer and Holmes (10). First, following Ref. (12), make a 
second change of variables 


v = x [S20] 

F u u = —F v x — uty, [S21] 

where the partial derivatives F u and F v are again evaluated at u = 0 and v = 0. This transformation converts the system of Eq. (S10) and 
Eq. (Sll) to 


x = -wy + f{x,y) [S22] 

y = tox + g(x, y) [S23] 

where 

f(x,y) = F[v(x),u{x)]+ujy [S24] 

g(x,y) = -^(F v F[v(x),u{x)] +F u G[v{x),u(x)]) -ojx. [S25] 
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102 Then compute 


103 a — ^ (/ill + fxyy + Qxxy + Qyyy) + [fxyifxx + fyy) [dxyidxx + flyi /) /xxQxx + 

104 Here the subscripts denote partial derivatives evaluated at (0,0). For example. 


[S26] 


105 


d 2 f 

dxdy 


( 0 , 0 ) 


[S27] 


106 When a < 0 the bifurcation is supercritical, and when a > 0 the bifurcation is subcritical (10-12). The calculation, which is performed using 

107 Mathematica, does not yield interpretable expressions. The results are instead expressed graphically by the color-coded stability plots in 
ioa Figures 3 and 6. 

109 Roughly speaking, instability results when the ballast feedback grows at a faster rate per mole than it is damped. The maximum growth rate 
no is of order y,8/c p whereas the damping constant is unity (in dimensionless time). The roughly hyperbolic stability boundary in the plane of 8 
in and y (Figure 3C) therefore corresponds to 8y ~ constant. Interpretation of the boundary with respect to the crossover concentration c x is 

112 somewhat more subtle. The per-mole growth rate of the ballast feedback is greatest when the COg~ concentration c = c x , where the slope of 

113 s(c, c x ) is steepest. Thus instability is most favored when c x is near c*. When c x is relatively far from c*, the growth rate is smaller and larger 

114 /i is required to outpace the homeostatic feedback. Consequently the boundary in Figure 3B is concave upwards, with its minimum at c x = c*. 


115 Size of the limit cycle. This section provides approximate results for the dependence of the amplitude of large limit cycles on the parameters 

116 of the model. The analysis follows from consideration of the large-c and low-c limits of the nullclines—the curves where corib vanish. Setting 

117 the right-hand side of Eq. (11) equal to zero, one finds the nullcline where c = 0: 

us w = wo — m[1 — bs(c, c p ) — 8s(c, c x ) — v\ = w cn ■ [S28] 

119 Likewise, Eq. (12) provides the nullcline where w = 0: 

120 w = wo + y [1 — bs(c,c p ) + 8s(c,c x ) + iy] = w w „. [S29] 

121 Size along the w-axis. The maximum value of w occurs when the trajectory crosses the low-c end of the w-nullcline. The maximum w can 

122 therefore be no greater than the low-c limit of equation Eq. (S29). Let R w to be the distance along the w-axis from this maximal excursion to 

123 the fixed point w*, normalized by fi: 

124 R w 

125 

126 The second relation used Eq. (15) and Eq. (S29). This quantity sets the scale of the maximum possible positive perturbation of w from the 

127 fixed point w* . To compare this to numerical simulation, define w max be the maximum observed value of w in a limit cycle. Then 

128 Aw = Wrnax — W* [S32] 

129 represents the maximum excursion relative to the fixed point for a particular choice of parameters. Figure SI 1 compares Aw//x, as computed 

130 for random choices of /x, 9, b , and c x , to R w . We see, as expected, that A w/ fi —> R w as R w increases. 


— (lim Wwn — w*') 
/x \c- >o J 


1 + 


eel 




[530] 

[531] 


131 Size along the c-axis. The magnitude of the change in w between the small-c limit of the w-nullcline and the large-c limit of the c-nullcline 

132 represents a characteristic scale of the change of w as c increases between crossings of the nullclines. Denoting it by R' w and normalizing by fi, 

133 one finds 


134 


135 


R' — - 

M 


lim Wen — lim w wn 

2 —Foo c —y 0 


= \2 + 9-b\. 


[533] 

[534] 


136 To obtain the characteristic change R c along the c-axis associated with R' w , write 


Rw I Vw | 

Rc Vc 


[S35] 


138 

139 

140 

141 

142 

143 

144 


where v m and v c are characteristic w- and c-components of velocity, respectively, in the phase plane. Although v w and v c can vary greatly 
during one period of a limit cycle, only their characteristic ratio as the trajectory travels between the maximum in w and the maximum in c is 
required. To a reasonable approximation, this ratio can be viewed as a typical value of the buffer factor /(c), which should be of the order of 
fo. Then 


K 1 

Rc ~ /o’ 


[S36] 


which yields 


Rc ~ /o|2 + 8 — b\. 


[S37] 
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145 The carbonate concentration c reaches its maximum c max when it crosses the c-nullcline w cn for large c. The quantity 

146 Ac = Cmax - C* [S38] 

147 is then the maximum excursion of c with respect to the fixed point. Figure S12 compares to Ac to R c . The nonlinear shape of the function 

148 likely indicates a limitation of the validity of Eq. (S36). The constant /o will always overestimate the characteristic value of /(c), but as limit 

149 cycles become larger and spend more time at high values of c, the approximation should improve. Thus the slope of the function steepens 
iso towards unity. 
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Fig. SI. The piecewise burial function (red dashed lines) of Zeebe and Westbroek (4) compared to the sigmoid function (blue curve) given 
by Eq. (7), with 7 = 4 and c p = 110 /rmol kg -1 . Parameters defining the piecewise function are taken from Table 1 of Ref. (4), under the 
assumption that neritic (shallow water) CaCC >3 production is 25% of the total production when [CO 3 - ] = 200 /xmol kg -1 . The total output is 
normalized to unity. 
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Fig. S2. Schematic diagram of the ballast-feedback model, a special case of the more general model illustrated by Figure 2. Fluxes are 
indicated in units of ji n . 
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Fig. S3. The partial derivatives dc/da and —dc/dw (dotted lines) as a function of c, computed for w = 2280 /imol kg \ 2.5 °C, 350 bar, 
and 35% salinity (1). The smooth curve is the approximation given by Eq. (TO) of the main text, for fo = 0.694, c/ = 43.9 /rmol kg -1 , and 
P = 1.70. 
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Fig. S4. Time series corresponding to the limit cycle of Figure 4 of the main text. A and B depict the variables w and c, respectively, obtained 
by integrating Eqs. (11-12) of the main text. The CO 2 concentration (C), the bicarbonate concentration (D), the pH (E), and the total alkalinity 
(F) are obtained from w and c by solving for the equilibrium of the carbonate system at each instant of time, assuming a temperature of 3.6 °C, 
a salinity of 35%, and a pressure of 350 bar (1). The dimensionless CO 2 feedback (G) and the dimensionless rate of carbonate burial (H) are 
computed from c(t). Time is given in units of the homeostat timescale r w . 
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Fig. S5. The dimensionless period T as a function of the rescaled amplitude A of the limit cycle, where A = (AibAc) 1 ^ 2 . The circles represent 
data obtained front numerical simulations using 531 random combinations of fi, 0, b, and c x that yield a limit cycle with Aw > 10 /xmol kg -1 . 
Here 2 < b < 8, 2 < 6 < 8, 150 < fi < 500 /.trnol kg -1 , and 65 < c x < 90 /rmol kg -1 . Fixed parameters are given in Table SI. The 
straight line is the prediction of Eq. (18) with Tq = 1. 
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Fig. S6. The maximum excursion of the DIC concentration w with respect to as a function of the dimensionless CO 2 injection rate v. As 
in Figure 5, c x = 55 fimo\ kg -1 and other parameters are given in Table SI. The critical injection rate v c ~ 0.391. 
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Fig. S7. Excitation near the saddle-node bifurcation of cycles. The phase space trajectory (A) and time series (B) are obtained for an 
above-threshold excitation with c x = 55.73 /xmol kg -1 . All other parameters are the same as in Figure 5C-D. Only the fixed point is stable. 
However, because the system is close to the saddle-node bifurcation of cycles, the excitation results in a series of pulses of similar magnitude 
before returning to the fixed point. 
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Fig. S8. The dimensionless size Aw //r (blue dots) of near-threshold excitations in the carbon cycle model as a function of the dimensionless 
time t/t w from the onset of an excitation to its peak. The red line is a rescaled graph of Aw = aji n r with a = 1.77. Simulations are 
performed with varying combinations of c x and 9 that lie in the excitation region of Figure 6 where v c < 1. Remaining parameters are 
specified in Table S1 . 
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Fig. S9. Dependence of the short-time threshold v' c on the injection time ti in numerical simulations of the carbon-cycle model. Inset: 
Comparison to the prediction of Eq. (24) on logarithmic axes. Parameters are the same used in Figures 5C and S6. 
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Fig. S10. Steady-state solutions. (A) The fixed point c* (Eq. (14)) as a function of c v and b. (B) The fixed point w* (Eq. (15)) as a function of 
c x and 6. Values of unspecified parameters are given in Table SI. 
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Fig. S11 . Comparison of Aw/^i, the /r-normalized difference between the maximum level of DIC during a limit cycle and the level of DIC of 
the unstable fixed point w*, to the prediction R w of equation Eq. (S31). The data points represent 531 random combinations of fi, 9 , b , and 
c x that yield ru ma x — w* > 10/rmol kg” 1 . Here 2 < b < 8, 2 < 9 < 8, 150 < /r < 500/rmol kg 1 , and 65 < c x < 90/rmol kg 1 . The 
identity line represents equality between the measurement and the prediction of equation Eq. (S31). The theory is asymptotically accurate for 
large limit cycles. (The Spearman rank correlation coefficient is 0.99.) 
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Fig. SI 2. Comparison of Ac/fi, the rescaled difference between the maximum level of COg - (c max ) during a limit cycle and the level of 
COg - of the unstable fixed point c*, to the scaling of equation Eq. (S37). The data points represent random combinations of the parameters as 
indicated in Figure S11. The Spearman rank correlation coefficient is 0.97. The curvature reflects changes in the characteristic value of the 
buffer factor /(c) as the size of the limit cycle increases. 
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Table SI. Parameters of the carbon-cycle model (Eqs. (11-12)) 
and their default values. 


symbol value units description 



250 a 

yimol kg -1 

b 

4 

— 

e 

5 

— 

c x 

70 

ytmol kg -1 

Cp 

110 b 

yimol kg -1 

IS 

0 

— 

Wo 

2000 

ytmol kg -1 

7 

4 b 

— 

fo 

0.694 c 

— 

c f 

43.9 C 

yimol kg -1 

p 

1.70° 

— 


a Obtained from the product of ji n — 0 
yr (7-9, 13). 

“From Figure S1. 

“From Figure S3. 


characteristic concentration ji n r w 
maximum CaC03 burial rate 
maximum respiration feedback rate 
crossover [CO^] (respiration) 
crossover [CO^j (burial) 

CO2 injection rate 
reference DIC concentration 
sigmoid sharpness index 
maximum buffer factor 
crossover [CO| — ] (buffering) 
sigmoid sharpness index 

.025 //mo! kg -1 yr -1 (6) and r„. — 10 4 


18 of 19 


Daniel H. Rothman 



151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 


References 


1. Zeebe R, Wolf-Gladrow D (2001) CO 2 in Seawater: Equilibrium, Kinetics, Isotopes. (Elsevier Science). 

2. Egleston ES, Sabine CL, Morel FMM (2010) Revelle revisited: Buffer factors that quantify the response of ocean chemistry to changes in 
DIC and alkalinity. Global Biogeochemical Cycles 24(1). GB1002, doi:10.1029/2008GB003407. 

3. Sarmiento JL, Gruber N (2006) Ocean Biogeochemical Dynamics. (Princeton, Princeton, N. J.). 

4. Zeebe RE, Westbroek P (2003) A simple model for the CaC03 saturation state of the ocean: The “Strangelove,” the “Neritan,” and the 
“Cretan” Ocean. Geochemistry, Geophysics, Geosystems 4( 12): 1104. 

5. BroeckerW (2014) Wally’s Carbon World. (Eldigio Press). 

6. Cai WJ, et al. (2008) A comparative overview of weathering intensity and HCO^ flux in the world’s major rivers with emphasis on the 
Changjiang, Huanghe, Zhujiang (Pearl) and Mississippi Rivers. Continental Shelf Research 28(12): 1538 - 1549. 

7. Archer D, Kheshgi H, Maier-Reimer E (1998) Dynamics of fossil fuel CO 2 neutralization by marine CaCOs. Global Biogeochemical 
Cycles 12(2):259-276. 

8. Archer D (2005) Fate of fossil fuel CO 2 in geologic time. Journal of Geophysical Research 110:C09S05. doi: 10.1029/2004JC002625. 

9. Ridgwell A, Hargreaves J (2007) Regulation of atmospheric CO 2 by deep-sea sediments in an Earth system model. Global Biogeochemical 
Cycles 21(2):GB2008. 

10. Guckenheimer J, Holmes P (1983) Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. (Springer-Verlag, New 
York). 

11. Strogatz S (1994) Nonlinear Dynamics and Chaos. (Addison-Wesley, New York). 

12. Izhikevich EM (2007) Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. (MIT Press, Cambridge). 

13. Sundquist E (1990) Influence of Deep-Sea Benthic Processes on Atmospheric CO 2 . Philosophical Transactions of the Royal Society of 
London A: Mathematical, Physical and Engineering Sciences 331(1616): 155—165. 


Daniel H. Rothman 


19 of 19 



